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Q ! Abstract 

o 

We have simulated lattice QCD with 2 flavours of massless staggered quarks. 
An irrelevant chiral 4-fermion interaction was added to the standard quark 



action to allow us to simulate at zero quark mass. Thermodynamics was 



studied on lattices of temporal extent 6. Clear evidence for a second order 



chiral transition was observed and the critical exponents Pmag, 5, v and j m ag 
were measured. These exponents did not agree with those expected by stan- 
dard universality arguments. They were, however, consistent with tricritical 
behaviour. The tt and a screening masses were measured and showed clear 

evidence for chiral symmetry restoration at this transition. 
12.38.Mh, 12.38.Gc, 11.15.Ha 
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I. INTRODUCTION 



In standard methods for simulating lattice QCD, the Dirac operator (matrix) 
becomes singular in the limit of zero quark mass. Since these methods require inversion of the 
Dirac operator they therefore fail for zero quark mass. The iterative methods for performing 
these inversions become progressively less efficient, requiring a number of iterations which 
diverges as this zero mass limit is approached. This means that, even at the mass of the 
physical u and d quarks, simulations on modest size lattices become prohibitively expensive 
in computer time, and most simulations are performed with unphysically large light-quark 
masses. 

QCD with 2 light quark flavours is expected to undergo a second order chiral transi- 
tion at finite temperature, only in the zero mass limit. It would thus be very helpful to be 
able to simulate in this limit in order to extract detailed information about this transition 
such as critical exponents and the equation of state in the neighbourhood of this transition. 
Here we note a second difficulty with the standard action. In the zero mass limit on a finite 
lattice, the chiral condensate which we need to use as an order parameter vanishes on a 
configuration by configuration basis, which restricts how small a mass we can use on a given 
lattice size. This contrasts with the behaviour of most lattice spin models, where the magne- 
tization (the order parameter for this class of models) has a finite value close to the infinite 
volume limit (for large enough lattices), even for zero magnetic field, for each configuration. 
The way in which the magnetization averages to zero over the ensemble of configurations, 
which it must for zero applied magnetic field, is that the direction of magnetization changes 
from configuration to configuration to enforce this. On a large enough lattice, a reasonable 
approximation to the spontaneous magnetization is obtained by averaging the magnitude of 
the magnetization (or its square) over the ensemble. 

To enable simulations at zero quark mass, we have modified the standard staggered 
quark action by the addition of an irrelevant 4-fermion term which preserves the symmetries 
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of the lattice action . Such a term will not affect the continuum limit of the theory. It does, 
however, serve to render the Dirac operator non-singular at zero quark mass. In the hadronic 
phase, where the chiral symmetry is spontaneously broken, this is easy to understand. The 
chiral condensate provides an effective mass to the quarks through the 4-fermion coupling. 
Even where the 4-fermion coupling is small, as it must be close to the continuum limit, this 
"constituent" quark mass is large enough to keep the number of iterations required to invert 
the Dirac operator manageable. With the standard staggered fermion action, the finite time 
extent of the lattice keeps the Dirac operator from becoming singular in the plasma phase . 
It does, however, become nearly singular for topologically non-trivial gauge configurations 
due to the lattice remnant of the Atiyah-Singer index theorem ||. For our modified action, 
the conservation of the axial quark-number current is explicitly broken by the 4-fermion 
interaction, thus avoiding this consequence of the index theorem. This keeps the number of 
iterations required to invert the Dirac operator manageable in the plasma phase. 

Although early attempts to extract critical exponents from staggered-quark lattice 
QCD appeared promising H, later attempts showed significant deviations from the 0(4) 
or 0(2) exponents expected from universality arguments. |5],|5| Attempts to fit scaling with 
(3 = Q/g 2 and quark mass m to the scaling functions of the 0(4) and 0(2) 3-dimensional spin 
models failed for lattices of temporal extent N t = 4 While the situation appeared 

better for N t > 6 (N t = 6 was somewhat ambiguous) insufficient "data" existed to make 
this conclusive. 

Our earlier simulations with our new action on N t = 4 lattices at m = indicated 
that the transition was probably first order [Q. Ongoing simulations on an Nt = 4 lattice 
with a smaller 4-fermion coupling indicate that the chiral transition is either a first order 
or a very narrow second order transition. The results of our simulations on 12 3 x 6 and 
18 3 x 6 lattices summarized in an earlier letter [Hi] and presented in detail in this paper, 



indicate that the N t = 6 transition is second order, but with the critical exponents of the 
3-dimensional tricritical point, rather that those of the 0(4) or 0(2) spin models. We were 
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able to arrive at this conclusion because we were able to study the scaling of the chiral 
condensate with (3 at m = and thus extract the critical exponent f3 mag . 

Section 2 describes our modified action. In section 3 we discuss the critical phenomena 
relevant to these simulations. We describe our simulations and results in section 4. In section 
5 we discuss our fits to critical scaling. Section 6 gives our conclusions. 



II. LATTICE QCD WITH CHIRAL 4-FERMION INTERACTIONS. 



We have modified our staggered quark lattice QCD action by the addition of an 
irrelevant chiral 4-fermion interaction. In the continuum the Euclidean Lagrangian density 
of this theory is 

C = \F^F, V + m)V - ^[(W 2 - (fysraV) 2 ]- (1) 

Lattice field theories incorporating fermions interacting both through gauge fields and 



through quartic self- interactions have been studied before — see for example jTT|, |L2 
131-15 — with promising preliminary results. 



Equation [l] can be rendered quadratic in the fermion fields by the standard trick 
of introducing (non-dynamical) auxiliary fields a and 7r in terms of which this Lagrangian 
density becomes 



£ = ^F^F^ + 4>{p+ a + *7r 75 T 3 + m)V + ^| (a 2 + tt 2 ) (2) 



The molecular dynamics Lagrangian for a particular staggered fermion lattice transcription 
of this theory in which T3 is identified with £5, the flavour equivalent of 75, is 

L = -PY$- \M^nUUUU)) + ^ JT^W - £ IiV /7 (a 2 + tt 2 ) 

□ 6 6 s - 8 

+\ E(^ 2 + el + m + pj 2 + TO + \ E(^ 2 + * 2 ) (3) 

l s 

where 



A =0 + m + — + im) (4) 
16 Y 

with i running over the 16 sites on the dual lattice neighbouring the site on the normal 
lattice, e = (— iy+y+ z+t an d fl the usual gauge- covariant "d-slash" for the staggered quarks. 
The factor in front of the pseudo-fermion kinetic term is appropriate for the hybrid 
molecular dynamics algorithm with "noisy" fermions, where Aip are chosen from a complex 
gaussian distribution with width 1. The "dots" represent derivatives with respect to molec- 
ular dynamics "time" as distinct from normal time. We note that 7 = 3/A 2 . Although the 
determinant of A does not appear to be real, it becomes so in the continuum limit. Without 
the gauge fields, this theory reverts to the one studied in [14]], with 3Nf flavours. 



The advantage of this choice of the chiral 4-fermion interaction is that it preserves 
the axial U(l) chiral symmetry of the normal staggered quark lattice QCD action generated 
by 75^5 at m = 0. This means that, when chiral symmetry is spontaneously broken, the 
pion associated with £575 will be a true Goldstone boson and will be massless at m — 0, 
even for finite lattice spacing. Under this exact chiral symmetry the fields transform as 

ij){n) -> e-^ e(n) -0(n) (5) 

er(ra) + Z7r(n) -> e^\a(n) + m[n)\ (6) 

from which we find that 

Aij)(n) -> e^ £ W^(n) (7) 

a{n) + te{n)n{n) -> e iMn) [a{n) + ie(n)ir(n)], (8) 

when m = 0. Hence, for massless quarks the above Lagrangian has an exact £7(1) axial 
flavour symmetry, the same remnant of continuum chiral flavour symmetry as the standard 
staggered lattice Lagrangian. 
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III. CRITICAL BEHAVIOUR OF THE CHIRAL TRANSITION IN LATTICE 

QCD 



Near a critical point, the order parameters of the theory in question scale as powers 
of the parameters in the theory. Since chiral symmetry is restored at the finite temperature 
transition from hadronic matter to a quark-gluon plasma, the appropriate order parameter 
is the chiral condensate, $ip). The equation of state describing the scaling of this order 



parameter close to the critical point is of the form [IB 



(W) = h^f QCD {th- l l^ s ) (9) 

where h = m/T and t — (3 — (3 C oc (T — T c )/T c . Jqcd is called the scaling function. This also 
serves to define 2 of the critical exponents, f3 mag and 5. Since {ipip) is non-zero for h = 0, 
t < equation | requires 

f Q c D (x) = C(-xf ; x^-oo, (10) 



so that 



(W) = C(-t) Pma '; t<0. (11) 



At t = 



(W) = Dh 1 / 5 ; D = f QCD (0). (12) 

Clearly, since (ipip) = 7(0"), similar scaling relations exist for (a). In addition, at h — 0, the 
a screening mass vanishes as 

m a = A\t\ v , (13) 

since the correlation length £ of the system diverges as £ oc \t\~ u . Finally we need to 
consider the scaling relations for the susceptibilities which measure the fluctuations in the 



order parameter (s). For reasons which will be explained later, we consider the scaling of the 
o susceptibility (Xa), rather than that for (ipip), with t. Near t — 0, Xc diverges as 

X, = V[((a) 2 )-((a)) 2 ]=c\tr^. (14) 

Here the inner angle brackets represent averages over the lattice for a single configuration 
while the outer angle brackets represent averages over the ensemble of configurations. 

Let us briefly give the arguments which indicate the expected universality class for 
this transition. First one notes that the Matsubara frequencies for bosonic excitations are 
2mrT, while those for fermionic excitations are (2n + 1)ttT [[17]]. The expectation that only 
low lying excitations contribute to critical phenomena leads one to the expectation that, 
at high T, only the lowest bosonic excitation n = contributes and the fermions decouple. 
Since this lowest bosonic mode is independent of time, this leads to an effective 3-dimensional 
theory (dimensional reduction) |18[ . Since the masses of the (tc, a) chiral multiplet vanish 
at the chiral transition, these excitations alone should determine the critical scaling at this 
transition. The effective field theories describing the interactions of the (tt, a) multiplet 
are 3-dimensional 0(4) sigma models. Thus 2-flavour QCD should lie in the universality 
class of the 0(4) sigma model in 3-dimensions with critical exponents f3 mag = 0.384(5), 
5 = 4.85(2), v = 0.749(4) and j mag = 1.471(6) @. Of course, lattice QCD with staggered 
quarks has only one Goldstone pion, and this symmetry is reduced to 0(2). Thus if the 
lattice is coarse enough that this symmetry breaking is relevant, 2-flavour QCD should lie 
in the universality class of the 0(2) sigma model in 3-dimensions, with critical exponents 
P mag = 0.35(1), S = 4.81(1), v = 0.679(3) and lmag = 1.328(6) |J. 

Of course, other possibilities exist. First the transition could be first order, which 
appears to be the case at JVj = 4 with our first choice of 4-fermi coupling JJ. Once the 
lattice spacing is large enough that lattice artifacts can affect the transition, more sym- 
metry breaking and symmetry conserving operators come into play giving a more complex 
transition. The simplest of these is the tri critical point. The sigma models which give the 
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universal critical points discussed in the previous paragraph lie in the universality class of 
the polynomial sigma models which are quartic in the fields. Tricritical behaviour occurs 
in sigma models which are sixth order in the fields ||20| . This means that there is a second 
symmetry preserving parameter or temperature and a second symmetry breaking parameter 
or magnetic field. Since 3 is the upper critical dimension for such sigma models, the critical 
exponents are given by a mean field analysis. Thus /3 mag = |, |, the S = 5, 3, v — ^, 7 = 1, |. 
The 2 values of f3 magi 5 and 7 describe the scalings when the 2 'magnetic fields' are varied 
independently. Even more values are possible when we take into account the second 'tem- 
perature'. A tricritical point occurs where a line of first order phase transitions continues 
as a line of second order phase transitions or in general where there are 2 competing phase 
transitions. 



IV. FINITE TEMPERATURE SIMULATIONS AND RESULTS FOR N T = 6. 



Our Nt — 4 simulations at 7 = 10 which showed evidence for a first order transition 
were reported in an earlier paper [0]. Here we present the details of our simulations on 



12 x 6 and 18 x 6 lattices. A brief summary of these results were given in a letter |10 
Our 12 3 x 6 simulations were performed at 7 = 10 and 7 = 20 while our 18 3 x 6 simulations 
were done at 7 = 20. 

We shall first discuss our results for zero mass quarks. Here, without the benefit of 
the symmetry breaking interaction (quark mass term), the direction of symmetry breaking 
changes over the run so that (a) and (it) and the corresponding chiral condensates (ipip) 
and (tp'js^ip) (where the () is averaging only over the lattice, not the ensemble) average 
to zero over the ensemble. It is in this way that the absence of spontaneous breaking of 
this continuous symmetry on a finite lattice is enforced. Note that this contrasts with the 
standard lattice formulation, where if we could simulate at zero quark mass, the chiral 
condensates would vanish on a configuration by configuration basis, and there would be no 
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way to predict the infinite lattice results. In our case, we can estimate the infinite lattice 
results by averaging the magnitudes yj {a} 2 + (n) 2 and \J '{ipip) 2 — ('^Ts^V') 2 , introducing an 
error which is of order 1/V below the phase transition and 1/y/V above it, on a lattice of 
4-volume V. 

The simulations of equation |^ were performed using the hybrid molecular-dynamics 
algorithm with "noisy" fermions, allowing us to tune the number of flavours to 2. The 
simulations on the 18 3 x 6 lattice were performed at 14 (3 values from 5.39 to 5.45. The 
behaviour of the chiral condensate and Wilson line (Polyakov loop) are given in figure [I]. The 




5.38 5.40 5.42 5.44 5.46 

FIG. 1. '(V'V')' an d Wilson line as functions of (3 on an 18 3 x 6 lattice at 7 = 20. 
values of the chiral condensate and (a) (actually the magnitudes defined above) as well as the 
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length of each run in molecular-dynamics time units are given in table |. (Note this "time" 
is twice that of equation |3| for consistency with that used by the HEMCGC and HTMCGC 
collaborations.) Let us make several observations with regard to these measurements. The 



P 


win 




'Time' 


5.3900 


0.4630(18) 


0.02543(10) 


5000 


5.3950 


0.4430(32) 


0.02454(16) 


5000 


5.4000 


0.4186(23) 


0.02337(11) 


10000 


5.4050 


0.3990(33) 


0.02244(17) 


10000 


5.4100 


0.3719(23) 


0.02109(12) 


20000 


5.4150 


0.3259(58) 


0.01880(30) 


40000 


5.4175 


0.2973(74) 


0.01734(40) 


50000 


0.4200 


0.2491(89 J 


0.014/9(48) 


50000 


5.4225 


0.1535(96) 


0.00957(53) 


50000 


5.4250 


0.0973(81) 


0.00649(45) 


20000 


5.4300 


0.0664(33) 


0.00479(18) 


20000 


5.4350 


0.0542(25) 


0.00414(13) 


10000 


5.4400 


0.0439(12) 


0.00364(7) 


10000 


5.4500 


0.0353(8) 


0.00326(5) 


5000 




TABLE I. (tpifj), (cr) and run lenj 


5th versus on an 18 3 x 6 lattice at 7 = 20. 





first is that the phase transition occurs at 5.4225 ^ j3 c ^ 5.4250. The second is that that 
due to fluctuations both (i/ji/j) and (a) are appreciable in the plasma phase. The third 
observation is that there are considerable deviations from the expected relation 

= (15) 

Part of this is attributable to these quantities having different finite volume errors, which is 
why the deviations become worse as (3 increases. The second is 0(dt 2 ) (and higher) errors 
which are quite large at this 7 for the dt = 0.05 we use. These errors affect (a) much 
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more than (ipip) and could presumably have been reduced by scaling the 'kinetic' term in 
equation |3] by a scale factor greater than 1, thus slowing the time evolution of a and ir. 
For this reason we perform our critical scaling analyses on (^ip) rather than (a). These 
assertions have been quantified for our N t = 4, 7 = 20 simulations and are presented in 
table |TJ, which suggests that both effects are coming into play. In fact, assuming a 1/V 
dependence for (?/>^>) and (a) would give ~ 0.7643 and ~ 0.03992 for their infinite volume 
limits on an N t = 4 lattice at /3 = 5.27 and dt = 0.05, giving a ratio of ~ 19.15. If, in 
addition, we assumed the coefficient of 1/V to be the same at dt = 0.0125 we would find 
('tpip) ~ 0.7567 and (a) « 0.03797 and the ratio ~ 19.93 at this dt. Assuming the finite dt 
error is 0(dt 2 ) leads to an estimate of ~ 19.98 for the infinite volume zero dt limit which 
is consistent with the correct value of 20. Our results for the 12 3 x 6 runs at 7 = 20 are 





8 3 x 4 lattice 






dt 






(^)/(a) 


0.0125 


0.7853(42) 


0.04051(22) 


19.38(4) 


0.0250 


0.7898(42) 


0.04114(21) 


19.19(4) 


0.0500 


0.7929(38) 


0.04246(19) 


18.68(4) 


0.1000 


0.8175(30) 


0.04715(16) 


17.34(3) 




12 2 x 24 x 4 lattice 






dt 






(W)/(a) 


0.0500 


0.7728(27) 


0.04067(14) 


19.00(3) 



TABLE II. dt and volume dependence of {i/)ip), (<r) and their ratios at = 5.27. 



presented in table |TlT] , while those for our 12 3 x 6 runs at 7 = 10 are given in table [IV]. 

The smoothness of the transition with no hints of metastability suggests a second 
order transition. Further evidence for this is the fluctuations with long time constants which 
are obvious from time histories of the chiral condensate close to the transition for the 18 3 x 6 
lattice at 7 = 20, shown in figures 0, |]. 
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& 




'Time' 


5.39000 


0.478(3) 


3500 


5.39500 


0.444(7) 


2000 


5.40000 


0.427(3) 


3500 


5.40500 


0.420(3) 


5000 


5.41000 


0.376(4) 


6500 


5.41500 


0.339(5) 


7500 


5.41750 


0.299(6) 


40000 


5.42000 


0.261(4) 


50000 


5.42125 


0.245(7) 


30000 


5.42250 


0.210(5) 


40000 


5.42500 


0.200(5) 


45000 


5.43000 


0.137(6) 


75000 



TABLE III. (t^ip) and run length versus (3 on an 12 3 x 6 lattice at 7 = 20. 



On the largest (18 3 x 6) lattice, we have measured the screening propagators for the 
7r and a fields, i.e. the n and a propagators for spatial separations. For this we use the n 
and a auxiliary fields. This has the advantage over using fermion bilinears with the same 
quantum numbers in that it automatically includes the disconnected contributions to both 
these fields. Below the transition, we define our a field by globally rotating our fields so that 
(71) = on each configuration. This is the best we can do since the direction of symmetry 
breaking rotates from configuration to configuration. Our zero momentum a and n screening 
propagators are then defined by 



v z 



^2<r(x,y,z,t) J2 v(x',y',z + Z,t') 

xyt x'y't' 

J2n(x,y,z,t) n(x',y',z + Z,t') 

xyt x'y't' 



- N x N y N t {af 



(16) 
(17) 



We note that P(Z), so defined, obeys 
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p 




'Time' 


5.4 


0.648(5) 


1500 


5.41 


0.616(4) 


1500 


5.42 


0.585(5) 


1500 


5.43 


0.539(2) 


10000 


5.44 


0.496(2) 


10000 


5.45 


0.436(2) 


35000 


5.455 


0.389(2) 


70000 


5.46 


0.330(3) 


45000 


5.46125 


0.302(3) 


60000 


5.4625 


0.267(3) 


60000 


5.46375 


0.234(4) 


35000 



(18) 



TABLE IV. (tpip) and run length versus j3 on an 12 3 x 6 lattice at 7 = 10. 

z 

Hence we should fit P(Z) to the form 

P(Z) — > Alexp^—Zm^) + exp(—(N z — Z)m- K )] + constant 



(19) 



Such fits in this low temperature domain yield values of m n in the range 0.1 < m v < 0.3. 
However, good fits are also obtained by fitting to a zero mass pion propagator. The reason 
for this is clear. Over the relatively small range we can fit these propagators on an N z = 18 



lattice, equation |19] is well approximated by the parabola 



P(Z) -> Aml(Z - N z /2f + [2Aexp(-m 7T N z /2) + constant] 



(20) 



which is the form of the zero mass propagator. Hence when we graph the 7r mass, we will 
use our knowledge that it must vanish when the chiral symmetry is broken to set it to zero 
below the transition. We fit S(Z) to the form 
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18 3 x6 LATTICE /3 = 5.42 y=20 
0.4 i 1 




Q Q I I I I I I I I I I I I I I I I I I I I I I I I I I 

10000 20000 30000 40000 50000 

TIME 

FIG. 2. Time history of (tptp) at (3 = 5.42. 
S(Z) — > B[exp(— Zm a ) + exp(—(N z — Z)m a )] + C[exp(— Zm n2 ) + exp(—(N z — Z)m n2 )] 

(21) 

and a similar fit with C — 0, and choose the fit with the best signal in this low-temperature 
phase. In the high temperature phase the restoration of chiral symmetry means that we can 
no longer distinguish a and ir. The best propagator to use in this plasma phase is PS(Z) 
defined as 

PS(Z) = P(Z) + S{Z) (22) 
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18 3 x6 LATTICE £ = 5.4225 7=20 
0.4 i 1 




Q Q I I I I I I I I I I I I I I I I I I I I I I I I I I 

10000 20000 30000 40000 50000 

TIME 

FIG. 3. Time history of (ijnf>) at p = 5.4225. 
Note here we do not have to rotate our fields. We then fit to the form 

PS{Z) -> A[exp(-Zm n ) + exp(-(N z - Z)m«)\. (23) 

The screening masses obtained in this way are plotted in figure |j. We note the a mass 
approaches zero as we approach the transition from the low temperature side. The combined 
a/ 7i mass rises from zero as we enter the plasma phase. However, even at the highest (3 
and hence T we consider, it is still far below the expected infinite temperature limit of 
2ttT = 7r/3. 

In addition, we measured contributions to the energy density, pressure and entropy 
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18 d x6 Lattice 7=20 



1.25 - 



1.00 



0.75 - 



0.50 - 



0.25 



0.00 




5.38 



5.40 



5.42 



5.44 



5.46 



P 



FIG. 4. a and 7r screening masses as functions of j3. Note that we have forced the the pion 
mass to be zero below the transition and the a and tt masses to be equal above the transition. The 
dashed line is at ir/3. 

density. Extraction of energy density or pressure would require measurements at the same 
(3 values at zero temperature. Hence we concentrate on the entropy density which can be 
determined from our measurements. To 1-loop order the partial energy densities for the 
gluons and the quarks are given by [511] 



^ = 4 4_6_ 

T 3 3 V 



i + sr 



1 (dC T dC a 

2 [~df~~dr. 



[(Pss) - (Pst)} 



(24) 
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and 



T 3 3*11^/4 



i / # 2 , 2 aa 



2 ' ' Nf. 2 2v 

7-^(a 2 + 7r 2 ) - 1 



2 V 2tt 2 3 «9£ 



(25) 



where P ss and P st are the space-space and space-time plaquettes respectively. The coefficient 
functions are given in reference . Since such expansions in terms of lattice couplings are 



notoriously poor, and an improved staggered fermion lattice perturbation expansion |22 
for these quantities has yet to be performed, we present graphs of the tree level as well 
as the 1-loop results for these quantities in figures |5], ||. The Stephan-Boltzmann value for 
s g (iue)/T 3 is 327r 2 /45 = 7.0184. The finite lattice size/spacing effects on an 18 3 x 6 lattice 
increase this to 8.0951. The Stephan-Boltzmann value for s/( erm j)/T 3 is 7ir 2 Nf/l5 = 9.2116 
which becomes 16.9823 on an 18 3 x 6 lattice. We note that both partial entropies increase 
rapidly from rather small values near the critical point and flatten out as the system passes 
into the plasma phase. The gauge field partial entropy lies above the Stephan-Boltzmann 
value in the plasma region, while the quark field partial entropy is about half its Stephan- 
Boltzmann limit. However, it should be noted that, since the njo screening masses indicate 
that the quarks in these mesonic excitations are still deeply bound at these (3 values, the 
high temperature limit would not be expected until a considerably higher (3. 

In addition to these zero mass simulations, we have performed simulations at finite 
mass on a 12 3 x 6 lattice. For these simulations we fixed j3 = (3 C . Here, the presence of even 
a relatively small symmetry breaking mass term was sufficient to align the condensate in the 
direction of (ipip) (or (a)) so that we were able to measure the chiral condensate directly. 
The chiral condensate at 7 = 10 as a function of quark mass is given in table 0, while that 
for 7 = 20 is given in table [VT1. 
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FIG. 5. Gluon partial energy density as a function of (3 
V. CRITICAL SCALING 



Much of what is presented here was presented in our letter JTIJ. We repeat these 
results here for completeness, along with results which have not been presented before. 

Fitting the measured values of ('ipip) on a 7 = 20, 18 3 x 6 lattice given in table | to the 
scaling form ($1/;) = C(J3 C - (3) Pma! > for the interval 5.41 < (3 < 5.4225 yields (3 C = 5.4230(2), 
C = 1.19(12) and j3 mag = 0.27(2) with a 94% confidence level. On a 12 3 x 6 lattice, also at 
7 = 20, fits to the values in table 011] gave (3 C = 5.4249(7), C = 1.46(22) and (3 mag = 0.32(4) 



18 



e 

Sh 

CD 



18 d x6 LATTICE y=20 
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FIG. 6. Quark partial energy density as a function of (5 

at a 27% confidence interval. Our 12 3 x 6 measurements at 7 = 10 (table [TV]) yielded 
(3 C = 5.4650(1) and (3 mag = 0.27(3). 

Clearly the sizable finite volume effects in these measurements are a major source 
of systematic errors. Having two different size lattices at 7 = enables us to remove the 
leading order finite volume effects. Here we note that 



«^) 2 - $7&ip) 2 ) = + constantly 



(26) 



where here the inner () stands for averaging over the lattice and the outer () stands for the 
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7 = 10, = 5.465 
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20000 
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TABLE V. Mass dependence of (tfiiff) at C for 7 = 10. 
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TABLE VI. Mass dependence of (if)i(>) at C for 7 = 20. 
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. Fitting 



this estimate of (ipip) to the scaling form over the range 5.41 < < 5.4225, used for the 
lattices of finite spatial volume gives C = 5.42270(15), C = 1.06(11), mag = 0.24(2) at 
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TABLE VII. Averages of squares of the chiral condensate on 12 3 x 6 and 18 3 x 6 lattices and 
the extrapolation to infinite spatial volume. 



a 78% confidence level. However, it is possible to fit this extrapolated chiral condensate 
over the complete range 5.39 < (3 < 5.4225, yielding f3 c = 5.42267(10), C = 1.00(3) and 
flmag = 0.229(9) while the confidence level only drops to 60%. The plot of this estimate of 
the chiral condensate with this latter fit superimposed is given in figure |7|. Finally we note 
that if one fixes f3 mag = 7, its tricritical value, one obtains a fit with confidence level 89% over 
the range 5.41 < (3 < 5.4225, which falls to 20% over the complete range 5.39 < (3 < 5.4225. 
From all these fits we conclude that the value of f3 mag is consistent with that of the 3- 
dimensional tricritical point, but inconsistent with either the 0(4) or 0(2) values. 
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FIG. 7. {4>ip) versus f3 with critical scaling curve superimposed. 

We now turn to consideration of the scaling of (ipip) with quark mass m at (3 = (3 C . 
This is difficult to perform with the usual staggered formulation, since it is hard to determine 
(3 C with sufficient precision. In addition, for m small enough, (ipip) vanishes linearly with 
m on a configuration by configuration basis at any (3. In our formulation, not only can (3 C 
be determined with adequate precision, but the chiral condensate remains finite on each 
configuration asm-»0 for (3 < f3 c and presumably should vanish as the appropriate power 
(1/5) of m at (3 = (3 C . What a finite mass does is orient the condensate in the required 
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direction. The value of m required to achieve such orientation should be somewhat smaller 
than that required to avoid the region of linear vanishing in the conventional approach. We 
fitted the 12 3 x 6, 7 = 10 results of table [V] to the scaling form (ijjip) = Dm 1 ^ . Our best 
fit gives D = 1.66(2) and 5 = 3.89(3). This fit has only a 4% confidence level. However, 
removing the point at m = 0.01 which lies 2.5 standard deviations off the curve increases 
the confidence level to 33% with negligible change to D and 8. We note that this 8 value is 
incompatible with the 0(4) or 0(2) values. It does, however, lie between the two 8 values 
(3 and 5) for the 3-dimensional tricritical point. The more general scaling can be obtained 
starting with the 6 effective field theory which describes tricritical behaviour. As argued 
in section 2, this can be treated in mean field approximation since it is at the upper critical 
dimension. The effective Hamiltonian for this theory is 

H = -d> 6 - mc-dj 5 + u-S 4 - mb-dj 3 + t-d> 2 - mad). (27) 
6^ 5^ V 3^ V V K ' 

The odd powers of dj are the symmetry breaking terms. Hence we identify them with the 
m(ijjijj) term in our Lagrangian, i.e. 

(^) = a0 + 6^ 3 + c^ 5 . (28) 

Note here that we have kept the third symmetry breaking operator <p 5 which has 8—1, 
since we do not have the luxury of redefining dj to remove it. Since we are considering the 
case where we are at the tricritical point for m = 0, the 2 reduced temperatures u and t are 
also zero. The vacuum expectation value of dj is that which minimizes H. Thus dj is the real 
positive solution of 

dJ> - m(a + bdj 2 + c 4 ) = 0. (29) 



Thus we fit (ipip) to the form of equation ^8| with dj given in equation ^9[ We obtained a 
fit with a = 1.05(2), b = 2.5(4), c = -3.1(1.4). Although this fit had only a 2% confidence 
level, removing the point at m = 0.01 gives a fit with parameters consistent with those 
including this point, but with a 35% confidence level. In fact, it is easy to convince oneself 
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that any reasonably smooth fit would have difficulty fitting this point. We have plotted 
the measurements of table [V] in figure |8] with this fit superimposed. Estimates of the errors 
induced in these measurements by the error in determining (3 C have been made from the 
tricritical scaling function, Jqcd and found to be small. 
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FIG. 8. Mass dependence of (iptp) at (3 C for 7 = 10. 



0.03 



Our measurements of the mass dependence of (ipip) on a 12 3 x 6 lattice at 7 = 20 
given in table [Vl| are somewhat poorer than those at 7 = 10. For this reason we are 
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unable to get reasonable fits to this 'data'. To clarify this, we have plotted this 'data' in 
figure |9|. Just looking at the points, we notice that the m = 0.004 and m = 0.005 points are 

7=20 
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FIG. 9. Mass dependence of (i>ip) at j3 c for 7 = 20. 

clearly incompatible. After several trial fits we found it that the m = 0.005 point should be 
excluded. Smoothness criteria caused us to also exclude the m = 0.02 point. Even then no 
fit came close to the m = 0.003 point. Here we suspect that this is because j3 c was more 
poorly determined than in the 7 = 10 case. Deviations of (3 from f3 c have the most effect on 
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the lowest mass point as is clear from equation |9[ Excluding these 3 points allows a fit with 
a = 0.83(3), b = 3.2(5), c = —2.5(1.5) at a 38% confidence level. This curve is plotted in 
figure |[ Clearly such a fit should be considered at best suggestive, so we base our conclusion 
that the scaling of (ipip) with mass is in accord with tricritical scaling purely on the 7 = 10 
measurements. 

From the plot of screening masses of figure it is clear that nothing quantitative 
can be obtained from the a masses in the low temperature domain. We have therefore 
concentrated our efforts on fitting the joint tt/<j screening masses in the plasma region to 
the scaling form m^/a = A(j3 — j3 c ) u . Since different treatment of these screening masses 
in the hadronic matter and quark-gluon plasma phases could bias our results, we chose to 



adopt the joint fit of the a and n propagators of equation |23] over the whole range of (3. This 
fit will be correct in the plasma phase. In the low temperature phase, the non-zero vacuum 
expectation value of the a field should drive the fitted 'mass' to zero. Fitting these screening 
masses to the screening form for (3 > /3 c (in fact we use all points for which (3 > 5.4225) we 
find A = 2.66(71), (3 C = 5.4196(7) and v = 0.59(7) admittedly with only a 0.3% confidence 



level. Plotting this fit on a graph figure 110 of mass measurements, it is clear why: the point 



at (3 = 5.35 is clearly inconsistent with those at (3 = 5.44 and (3 = 5.45. Removing the 
point at (3 = 5.35 improves the confidence level to an acceptable 39% while only changing 
the fitting parameters within errors (for example v = 0.62(6)). We attribute this to the 
systematic errors in performing fits to point-source propagators with a lattice extent of only 
18 and limited statistics. From these fits we conclude that v is consistent with the tricritical 
value (|), but does not rule out an 0(4) or 0(2) value. 

We have extracted a susceptibility from our o 'data' on the 18 3 x 6 lattice at 7 = 20 
using x<t — ^[((a) 2 ) ~~ {( cr )} 2 ]- Since our estimate of ((a)) was obtained after rotating 
(<t) to where (tt) = configuration by configuration this actually defines a 'radial' rather 
than the conventional susceptibility. However, this should diverge in the same manner as the 



conventional susceptibility, which is all that matters. Fitting to the form \o = c\f3 — f3 c 



'Jmag 
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FIG. 10. Joint a/ir screening masses as functions of (3. 

over the range 5.405 < (3 < 5.45 gives (3 C = 5.4215(1), c = 0.0096(11) and -f mag = 0.66(3) 
with a 31% confidence level. Our measurements and fit are plotted in figure |TT|. Since 
the value of this critical exponent lies between the 2 7 ma9 values (| and 1) associated with 
the tricritical point, but far from the 0(4) and 0(2) exponents this is further evidence for 
tricriticality. The reason we did not use the (ipip) susceptibility is that our measurements 
of (ipip) are noisy estimators using a single noise vector per configuration which gives extra 
contributions to the susceptibility estimate that can only be removed by use of more than 
one noise vector per configuration. 
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FIG. 11. a susceptibility function of (3. 

VI. DISCUSSION AND CONCLUSIONS 

We have presented evidence that the N t = 6 finite temperature phase transition for 
2-flavour staggered lattice QCD with extra chiral 4-fermion interactions is in the universality 
class of a 3- dimensional tricritical point rather than that of the 3-dimensional 0(4) or 0(2) 
sigma model, which was expected. Since a tricritical point occurs when there are extra 
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relevant operators with consequent competing transitions, we believe this to be due to the 
extra higher- dimensional operators which can be important on such coarse lattices. For this 
reason, we believe the 0(4)/ 0(2) universality class would be observed for sufficiently fine 
lattices, i.e. those with sufficiently large N t . Since the N t = 4 transition with 7 = 10 was 
observed to be first order, and a tricritical point marks the transition from first order to 
universal second order behaviour, this strongly suggests that 0(4)/ 0(2) universality would 
be observed for Nt as small as 8. 

The fact that our modified action allowed us to work at zero quark mass was crucial 
in allowing us to come to these conclusions. This enabled us to measure the critical exponent 
Pmag directly. Our best estimates lay in the range 0.23 < f3 mag < 0.27 (with errors typically 
around 2 in the second decimal place) which is consistent with the \ of the tricritical point, 
but inconsistent with the 0.384(5) of the 0(4) sigma model, the 0.35(1) of the 0(2) sigma 
model and the ~ of mean field theory. Because these scaling fits give precise estimates for 
P c , it is possible to measure 5 directly from the scaling of the chiral condensate with quark 
mass with (3 fixed at (3 C . This value, 5 = 3.89(3), lay between the two 5 values (3 and 5) 
of the tricritical point, and a satisfactory fit to the more complex scaling this implies was 
obtained. Again this behaviour was incompatible with 0(4)/0(2) scaling where S ~ 4.8 and 
mean field scaling with 5 = 3. The scaling of the a/n mass in the plasma domain yielded a 
v value which was compatible with tricritical, 0(4)/ 0(2) and mean field scaling. Finally we 
were able to extract the susceptibility critical exponent 7 mag = 0.66(3) which lies between 
the two values (| and 1) of the tricritical point but disagrees with the 0(4)/0(2) values 
(1.328(6) and 1.471(6)) and the mean field value (1). 

Since we believe we are seeing more complex phase structure which is driven by lattice 
artifacts, the behaviour we are seeing does not have to apply to the standard staggered 
action. However, since scaling analyses with the standard action on N t = 4 lattices showed 
definite departures from 0(4)/ 0(2) universality @J§J§,§] it strongly suggests that some such 
behaviour is present. In fact, based on our earlier, N t = 4 results one of the collaborations 
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suggested || that these anomalies might be associated with the first order transition we 
reported in [Q]. Unfortunately there are not yet sufficient high precision measurements 
with the standard action at N t = 6 to tell if it too displays departures from 0(4)/ 0(2) 
universality 

Clearly, simulations at N t = 8 should be next on our agenda. We would also suggest 
that the addition of such 4-fermion terms should be considered for other actions, especially 
improved staggered actions. 
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